function LE = LEs(eps_val, params, X0, T_trans, T_le, dt_save, RelTol, AbsTol, numLE)
% computeLE 计算给定参数下系统的前 numLE 个 Lyapunov 指数 (例如 LE1, LE2)
% 在计算对数放大因子时做截断，避免出现 -Inf 的情况
%
% 输入参数 (略) ...
% 输出:
%   LE: numLE×1 向量, LE(1) 为最大 Lyapunov 指数, LE(2) 为次大 Lyapunov 指数

    if isrow(X0), X0 = X0'; end
    odefun = @(t, X) hrfnhr(t, X, eps_val, params);
    opts = odeset('RelTol', RelTol, 'AbsTol', AbsTol);

    % 1) 去除瞬态
    [~, X_tr] = ode45(odefun, [0, T_trans], X0, opts);
    X_main = X_tr(end,:)';

    % 2) 初始化扰动向量
    delta0 = 1e-8;
    n = length(X_main);
    Q = zeros(n, numLE);
    for i = 1:numLE
        v = randn(n,1);
        v = v / norm(v) * delta0;
        Q(:, i) = v;
    end
    % 先做一次正交化
    [Q, ~] = qr(Q, 0);

    sum_log = zeros(numLE,1);
    t_current = 0;

    while t_current < T_le
        t_end = t_current + dt_save;

        % 积分主轨迹
        [~, X_sol] = ode45(odefun, [t_current, t_end], X_main, opts);
        X_main_new = X_sol(end,:)';

        % 对每个扰动方向做积分
        new_vecs = zeros(n, numLE);
        for j = 1:numLE
            X_pert0 = X_main + Q(:, j);  % 初始状态 + 扰动
            [~, Xp_sol] = ode45(odefun, [t_current, t_end], X_pert0, opts);
            X_pert_new = Xp_sol(end,:)';
            diff_vec = X_pert_new - X_main_new;
            new_vecs(:, j) = diff_vec;
        end

        % 对 new_vecs 做 QR 分解
        [Qn, R] = qr(new_vecs, 0);

        for j = 1:numLE
            % 原始放大因子: abs(R(j,j)) / delta0
            amp = abs(R(j,j)) / 1e-8;  % 因为之前扰动大小是 delta0=1e-8
            if amp < 1e-15
                amp = 1e-15;  % 避免 log(0) => -Inf
            end
            sum_log(j) = sum_log(j) + log(amp);
        end

        % 更新 Q
        Q = Qn;
        % 归一化: 确保新的 Q 列向量大小回到 delta0
        for j = 1:numLE
            normQj = norm(Q(:, j));
            if normQj < 1e-15
                % 若真的出现极小值，重置为随机向量
                rvec = randn(n,1);
                rvec = rvec / norm(rvec) * delta0;
                Q(:, j) = rvec;
            else
                Q(:, j) = Q(:, j) / normQj * delta0;
            end
        end

        % 更新主轨迹
        X_main = X_main_new;
        t_current = t_end;
    end

    % 计算 LE
    LE = sum_log / T_le;
end
